Charles Vin (21216136) & Aymeric Delefosse (21213744) DAC 2023-2024

TME 10- LIME / SHAP

L'objet de ce TME est de tester les algorithmes LIME et SHAP sur des tâches de classification assez simple.

import lime.lime_tabular
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import shap
from lime.lime_tabular import LimeTabularExplainer
from sklearn import datasets
from sklearn.datasets import make_moons
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import DecisionBoundaryDisplay
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

1 - Préparation du dataset

Utiliser le dataset UCI ML Breast Cancer Wisconsin (Diagnostic), comme dans le TME 3 et apprendre dessus un modèle linéaire de régression logistique

breast_cancer = datasets.load_breast_cancer()

X_train, X_test, y_train, y_test = train_test_split(
    breast_cancer.data, breast_cancer.target, test_size=0.25, random_state=42
)
feature_names = breast_cancer.feature_names
target_names = breast_cancer.target_names

logreg = LogisticRegression(max_iter=2500)
logreg.fit(X_train, y_train)
y_pred = logreg.predict(X_test)
print("=" * 50)
print(f"Accuracy = {(y_pred == y_test).sum() / len(y_test):.4f}")
print("=" * 50)

==================================================
Accuracy = 0.9650
==================================================

Samples per class: 212(M),357(B)

(breast_cancer.target == 0).sum(), (breast_cancer.target == 1).sum()

(212, 357)

breast_cancer.target_names

array(['malignant', 'benign'], dtype='<U9')

2- Explication avec SHAP

Nous allons étudier une explication donnée par SHAP.

Q2.1 - Exécuter le bloc ci-dessous pour obtenir les explications globales de SHAP. Commentez.

explainer = shap.LinearExplainer(logreg, X_train, feature_names=feature_names)
shap_values = explainer.shap_values(X_test)

shap.summary_plot(shap_values, X_test, plot_type="bar", feature_names=feature_names)

Commentaire :

Les valeurs sur l'axe des abscisses représentent la valeur absolue moyenne de SHAP pour chaque caractéristique dans toutes les instances de l'ensemble de test. Plus la valeur est élevée, plus l'impact de cette caractéristique sur les résultats du modèle est important. Le diagramme à barres nous permet donc de déduire que :

• mean perimeter a l'impact moyen le plus élevé sur les résultats du modèle. Cela suggère que le périmètre moyen d'une zone tumorale a une grande influence sur les prédictions du modèle, ce qui est probablement en corrélation avec la classification des tumeurs malignes ou bénignes.
• Les caractéristiques mean area et worst area ont également un impact significatif, ce qui est logique puisqu'elles sont liées à la taille de la tumeur, qui sont essentielles pour diagnostiquer le cancer du sein.
• Les caractéristiques telles que area error, mean radius, worst texture, mean texture, texture error et permimeter error ont un impact modéré.
• Les caractéristiques ayant le moins d'impact sur les résultats du modèle sont worst radius, worst concavity et worst compactness (et toutes les caractéristiques suivantes).

Cependant, ce graphique de synthèse ne nous indique pas la direction de l'impact (si la caractéristique augmente ou diminue la probabilité que le modèle prédise une classe), mais plutôt l'ampleur des contributions des caractéristiques à la sortie du modèle. Pour cela, il faut utiliser un autre type de diagramme SHAP, comme le diagramme beeswarm ou le diagramme de décision, qui montre les valeurs SHAP individuelles pour chaque instance.

De plus, une deuxième remarque est, à quel point la corrélation entre plusieurs caractéristiques rend la valeur de SHAP similaire ? On peut effectivement penser que toutes les caractéristiques liés à la taille (perimeter, area) sont corrélées entre elles. Les valeurs de Shapley ne supposent pas l'indépendance. Mais, à cause du coût de calcul lié, il est courant de supposer l'indépendance, même si c'est en pratique très peu le cas (approximation Kernel SHAP).

Il est également important de noter que les caractéristiques ayant un faible impact sur le modèle peuvent néanmoins être importantes ; elles peuvent simplement avoir un rôle plus nuancé ou interagir avec d'autres caractéristiques d'une manière qui n'est pas prise en compte par les seules valeurs moyennes du SHAP.

Regardons la valeur moyenne non absolue de SHAP pour chaque caractéristique.

shap_mean_df = pd.DataFrame(
    {"Feature": feature_names, "mean(SHAP value)": shap_values.mean(0)}
)

# Sorting the DataFrame by the SHAP values
shap_mean_df = shap_mean_df.sort_values(
    by="mean(SHAP value)", ascending=False, key=lambda x: abs(x)
).head(15)

# Plotting the graph using seaborn
plt.figure(figsize=(10, 8))
sns.barplot(
    x="mean(SHAP value)", y="Feature", hue="Feature", data=shap_mean_df, palette="bone"
)
plt.title("Mean SHAP values for model features")
plt.xlabel("mean(SHAP value) (average impact on model output magnitude)")
plt.ylabel("Feature")
shap_plot = plt.gca()

shap_mean_df

	Feature	mean(SHAP value)
3	mean area	-1.022887
2	mean perimeter	0.886764
13	area error	0.370578
0	mean radius	-0.367687
21	worst texture	-0.347981
1	mean texture	0.097638
23	worst area	0.062468
12	perimeter error	-0.054158
22	worst perimeter	0.028523
11	texture error	0.022610
25	worst compactness	-0.011328
26	worst concavity	0.005120
28	worst symmetry	-0.004997
20	worst radius	-0.003722
6	mean concavity	0.002436

Commentaire :

• Les caractéristiques mean perimeter et mean area restent celles dont la valeur moyenne de SHAP est la plus élevée, ce qui indique leur forte influence sur les prédictions du modèle.
  • Une valeur plus élevée pour mean area augmente fortement la probabilité que le modèle prédise une tumeur maligne (classe 0).
  • Une valeur plus élevée pour mean perimeter augmente fortement la probabilité que le modèle prédise une tumeur bénigne (classe 1).
• La worst area n'apparait pas aussi haut que précédemment : de 6 en mean(|SHAP value|), on passe à 0.06 en mean(SHAP value), ce qui suggère qu'il y a autant de valeurs très positives ou très négatives en fonction des données (ce qu'on verra ensuite avec le beeswarm).

Ainsi, fournir les valeurs non absolues moyennes de SHAP permet de se faire une autre idée que le graphique précédent, mais n'est pas parfaite, au vu de comment SHAP fonctionne.

shap_values_ = explainer(X_test)
shap.plots.beeswarm(shap_values_)

Commentaire :

Ou, toujours en utilisant SHAP, on peut obtenir un graphe similaire mais beaucoup plus joli (le fameux beeswarm). La caractéristique distinctive de ce graphique est l'utilisation de points individuels pour représenter les valeurs spécifiques de chaque caractéristique pour les différentes prédictions. Chaque point correspond à une instance de la caractéristique dans le jeu de données, rendant ce graphique très informatif.

• Pour mean perimeter, worst area, area error et worst texture, des valeurs élevées contribuent significativement à augmenter la probabilité que le modèle prédise la présence de cancer, alors que des valeurs faibles tendent à prédire son absence.
• A l'inverse, pour mean area et mean radius, des valeurs élevées sont plutôt associées à l'absence de cancer, et des valeurs plus faibles à sa présence.
• Les caractéristiques listées après worst perimeter semblent avoir un impact minime sur les prédictions du modèle.

Q2.2 - Exécuter le bloc ci-dessous pour obtenir une explication locale avec SHAP du premier exemple de la base de test. Commentez.

# shap.initjs()
shap.force_plot(
    explainer.expected_value,
    shap_values[0, :],
    X_test[0, :],
    feature_names=feature_names,
    matplotlib=True,
)

Commentaire :

Le graphique SHAP "force plot" est une représentation visuelle de l'impact des caractéristiques sur la prédiction d'un modèle pour une instance spécifique. L'axe des abscisses représente une transformation log-odds de la probabilité prévue.

Voici un résumé des nos observations :

• $f(x)$ : Ce n'est pas la probabilité prédite $ p $ mais le logit de la probabilité, c'est-à-dire $ \log \left(\frac{p}{1-p}\right) $. Dans le cas de notre premier exemple de la base de test, la valeur de $ f(x) $ est positive et étant donné que $ \textrm{logit}^{-1}(\alpha) = \frac{1}{1 + exp(-\alpha)} \Leftrightarrow p = 0.86$ pour $\alpha = 1.81$, cela suggère une certaine certitude concernant la classification de notre patient.

• Mesure d'incertitude : Ainsi, plus la valeur de $f(x)$ est éloignée de zéro, plus la prédiction du modèle est certaine. Une valeur fortement négative indiquerait une forte conviction que la tumeur est "maligne", tandis qu'une valeur fortement positive indiquerait une forte conviction que la tumeur est "bénigne".

• Caractéristiques positives (rose) : Les caractéristiques worst texture, worst perimeter, area error, worst area, et mean perimeter ont un effet positif sur la prédiction du modèle, indiquant que les valeurs élevées pour ces caractéristiques sont associées à une prédiction de la classe "bénigne".

• Caractéristiques négatives (bleu) : Les caractéristiques mean area et mean radius ont un effet négatif sur la prédiction du modèle, suggérant que les valeurs élevées pour ces caractéristiques sont associées à une prédiction de la classe "maligne".

• Valeurs des caractéristiques : Les nombres sur les barres indiquent les valeurs réelles des caractéristiques pour l'instance spécifique analysée. Ces valeurs sont utilisées pour calculer l'impact sur le logit de la probabilité.

Lors du TME, nous nous demandions qu'est-ce que la base value (qui correspond à l'expected value) et à quoi correspond-elle. La base value correspond à une valeur de prédiction d'un modèle sans caractéristiques (feature-less model), c'est-à-dire l'espérance de notre modèle, $\mathbb{E}[f(x)]$.

Q2.3 - Comparer avec l'exemple suivant de la base de test.

shap.force_plot(
    explainer.expected_value,
    shap_values[1, :],
    X_test[1, :],
    feature_names=feature_names,
    matplotlib=True,
)

Commentaire :

Dans ce cas, la sortie du modèle est beaucoup plus négative (-17.35), ce qui suggère une forte confiance que l'instance appartient à la classe avec la probabilité $ p $ proche de 0 ("bénigne") ($ p = 0.999999970826353 $). Ce qui est particulièrement intéressant ici, c'est la magnitude extrême de $ f(x) $ comparée à l'autre exemple. Cela souligne l'importance de comprendre comment les valeurs individuelles de chaque caractéristique contribuent à la prédiction finale dans le contexte spécifique des données sur lesquelles le modèle a été formé.

Comparons cela avec l'analyse précédente :

• Les caractéristiques qui ont le plus d'impact sur la prédiction sont toujours représentées par la longueur des barres, mais cette fois, elles poussent la prédiction dans la direction opposée.
• mean radius, et mean area ont des contributions positives, ce qui signifie que ces trois caractéristiques "poussent" le modèle vers une prédiction "bénigne".
• À l'inverse, mean perimeter, worst area, et area error poussent la prédiction dans une direction négative ("maligne") et contrebalancent les effets des autres caractéristiques.

Malgré les contributions positives des caractéristiques roses, la prédiction globale du modèle penche fortement vers la classe "maligne". Cet exemple illustre la complexité de l'interprétation des modèles prédictifs, en particulier lorsqu'on utilise des outils d'explication tels que SHAP. Il est essentiel de bien comprendre la direction et la magnitude des contributions de chaque caractéristique pour interpréter correctement ces graphiques. Malgré la "beauté" de ce type de graphique, il y a besoin d'un temps d'apprentissage pour bien les lire...

3 - Expérimentations

• Comparer les résultats fournis par LIME et SHAP sur des données identiques et commenter les résultats
• Examiner les résultats par exemple sur les données half-moons dans un cas où une explication linéaire est inappropriée
• Générer d'autres explications globales avec SHAP
• Modifier votre propre implémentation de LIME pour utiliser le noyau SHAP et comparer au résultat fourni par SHAP [facultatif]

explainer_lime = LimeTabularExplainer(
    X_test,
    feature_names=feature_names,
    class_names=target_names,
    categorical_features=[],
    mode="classification",
)

exp = explainer_lime.explain_instance(X_test[0], logreg.predict_proba, num_features=10)
exp.as_pyplot_figure()
plt.show()

Commentaire LIME :

• Les barres rouges indiquent les caractéristiques qui contribuent à une classification de la classe "maligne" et les barres vertes indiquent celles qui contribuent à une classification de la classe "bénigne".
• La largeur des barres reflète l'importance de la contribution de la caractéristique à la prédiction de la classe. Par exemple, une grande barre verte pour worst area indique que, pour cette instance, lorsque la valeur de worst area se situe entre 490.65 et 684.50, cela diminue fortement la probabilité d'une prédiction de cancer (suggérant que c'est bénin).

shap.plots.bar(shap_values_[0], max_display=10)

Commentaire SHAP :

• Les barres bleues indiquent les caractéristiques qui contribuent négativement à la prédiction (vers la classe "maligne"), et les barres rouges indiquent celles qui contribuent positivement (vers la classe "bénigne").
• La longueur des barres et la valeur SHAP associée montrent la force de cette contribution. Par exemple, mean area a une contribution négative importante, ce qui signifie que, pour cette instance, une grande valeur de mean area pousse la prédiction vers la classe 0 ("maligne").

Comparaison LIME vs SHAP :

Il y a une cohérence globale dans l'explication des caractéristiques par LIME et SHAP :

• worst area, mean perimeter, worst texture, area error et worst perimeter contribuent toutes à classifier la donnée vers la classe positive. La différence est à quel degré elles le font, en fonction de l'explication donnée par SHAP ou par LIME. Par exemple, worst area a une importance plus grande dans LIME que dans SHAP.
• mean area et mean radius sont toutes les deux des caractéristiques communes de nos deux explications, et contribuent à classifier la donnée vers la classe négative. Ici, contrairement aux caractéristiques précédentes, LIME et SHAP attribuent un poids similaire.
• Il y a des caractéristiques dans l'explication de LIME qui ne sont pas dans l'explication de SHAP, et vice-versa.

La différence clé est que SHAP prend en compte les interactions entre les caractéristiques et attribue des valeurs SHAP basées sur une répartition équitable des contributions, alors que LIME génère une explication locale en approximant le modèle par un modèle linéaire simplifié.

Corrélation et SHAP

(Sur le dataset Breast Cancer)

X = pd.DataFrame(breast_cancer.data, columns=breast_cancer.feature_names)
y = pd.Series(breast_cancer.target, name="target")

# Concaténer les features et la cible pour calculer la corrélation
df = pd.concat([X, y], axis=1)

# Calculer la matrice de corrélation
corr_matrix = df.corr()

# Visualiser la matrice de corrélation avec une heatmap
plt.figure(figsize=(20, 20))
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="coolwarm", center=0)
plt.title("Matrice de corrélation des caractéristiques de Breast Cancer")
plt.show()

En fait, je ne suis pas sûr de comment bien montrer une potentielle corrélation entre les features corrélées et les valeurs de SHAP. :/

Update : En faisant la question "Générer d'autres explications globales", j'ai compris qu'on pouvait utiliser le scatter plot de SHAP pour montrer une corrélation ! J'explique plus en détail ce graphique dans la question "Générer d'autres explications globales".

On se rappelle que les caractéristiques ayant le plus d'impact sur notre modèle sont mean perimeter, mean area, worst area, area error, mean radius et worst texture. Ici, on s'intéressera uniquement aux caractéristiques qui font écho à la taille : on peut voir qu'elles sont toutes (très) fortement corrélées les unes entre les autres grâce à la matrice de confusion ci-dessus.

shap.plots.scatter(shap_values_[:, "mean perimeter"])
shap.plots.scatter(
    shap_values_[:, "mean perimeter"], color=shap_values_[:, "mean area"]
)
shap.plots.scatter(
    shap_values_[:, "mean perimeter"], color=shap_values_[:, "worst area"]
)
shap.plots.scatter(
    shap_values_[:, "mean perimeter"], color=shap_values_[:, "area error"]
)

Et les graphiques ci-dessus nous le confirment. Autre le fait que mean perimeter est un effet linéaire sur le modèle, on peut en conclure, grâce aux couleurs, que plus une aire est grande, plus la prédiction penchera en faveur de "bénigne", et vice-versa. Cela fait bien écho aux corrélations positives élevées entre ces caractéristiques.

Données non linéaires

def plot_boundaries(X, y, clf, obs=None, ax=None):
    """Plot the data and the decision boundary resulting from a classifier."""
    if ax is None:
        fig, ax = plt.subplots()
    x_min, x_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
    y_min, y_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
    DecisionBoundaryDisplay.from_estimator(clf, X, ax=ax, eps=0.5, cmap=plt.cm.RdBu)
    # Plot the training points
    ax.scatter(X[:, 0], X[:, 1], c=y, edgecolors="k")
    if obs is not None:
        ax.scatter(
            obs[0],
            obs[1],
            c="red",
            edgecolors="k",
            marker="o",
            label="Instance à expliquer",
            s=100,
        )
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)

    return ax

# Générer des données half-moons
X, y = make_moons(n_samples=500, noise=0.20, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

# Entraîner un classificateur
tree = RandomForestClassifier(n_estimators=100, random_state=42)
tree.fit(X_train, y_train)
y_pred = tree.predict(X_test)
print("=" * 50)
print(f"Accuracy = {(y_pred == y_test).sum() / len(y_test):.4f}")
print("=" * 50)

plot_boundaries(X, y, tree)
plt.title("Frontières de décision sur half moons avec un arbre de décision")
plt.show()

==================================================
Accuracy = 0.9760
==================================================

# SHAP
shap_explainer = shap.TreeExplainer(tree, feature_names=["feature 1", "feature 2"])
shap_values = shap_explainer.shap_values(X_test)

# LIME
lime_explainer = lime.lime_tabular.LimeTabularExplainer(
    X_train,
    feature_names=["feature 1", "feature 2"],
    class_names=["class 0", "class 1"],
    discretize_continuous=True,
)

# Sélectionner une instance à expliquer
i = 10
ax = plot_boundaries(X, y, tree, X_test[i])
plt.legend()
plt.title("Frontières de décision sur half moons avec un arbre de décision")
plt.show()

print(f"Predicted class: {y_pred[i]}")
print(f"Real class: {y_test[i]}")

shap_explanation = shap_explainer.shap_values(X_test[i])
lime_explanation = lime_explainer.explain_instance(
    X_test[i], tree.predict_proba, num_features=2
)

# Afficher les résultats SHAP
shap.force_plot(
    shap_explainer.expected_value[1],
    shap_explanation[1],
    X_test[i],
    feature_names=["feature 1", "feature 2"],
    matplotlib=True,
)
# Afficher les résultats LIME
lime_explanation.show_in_notebook()

Predicted class: 1
Real class: 1

# Sélectionner une instance à expliquer
i = 42

ax = plot_boundaries(X, y, tree, X_test[i])
plt.legend()
plt.title("Frontières de décision sur half moons avec un arbre de décision")
plt.show()

print(f"Predicted class: {y_pred[i]}")
print(f"Real class: {y_test[i]}")

shap_explanation = shap_explainer.shap_values(X_test[i])
lime_explanation = lime_explainer.explain_instance(
    X_test[i], tree.predict_proba, num_features=2
)

# Afficher les résultats SHAP
shap.force_plot(
    shap_explainer.expected_value[1],
    shap_explanation[1],
    X_test[i],
    feature_names=["feature 1", "feature 2"],
    matplotlib=True,
)
# Afficher les résultats LIME
lime_explanation.show_in_notebook()

Predicted class: 0
Real class: 0

# Sélectionner une instance à expliquer
obs = np.array([0, 0.54])  # frontière

ax = plot_boundaries(X, y, tree, obs)
plt.legend()
plt.title("Frontières de décision sur half moons avec un arbre de décision")
plt.show()


shap_explanation = shap_explainer.shap_values(obs)
lime_explanation = lime_explainer.explain_instance(
    obs, tree.predict_proba, num_features=2
)

# Afficher les résultats SHAP
shap.force_plot(
    shap_explainer.expected_value[1],
    shap_explanation[1],
    obs,
    feature_names=["feature 1", "feature 2"],
    matplotlib=True,
)
# Afficher les résultats LIME
lime_explanation.show_in_notebook()

divide by zero encountered in scalar divide

obs = np.array([1, -0.15])  # frontière

ax = plot_boundaries(X, y, tree, obs)
plt.legend()
plt.title("Frontières de décision sur half moons avec un arbre de décision")
plt.show()

shap_explanation = shap_explainer.shap_values(obs)
lime_explanation = lime_explainer.explain_instance(
    obs, tree.predict_proba, num_features=2
)

# Afficher les résultats SHAP
shap.force_plot(
    shap_explainer.expected_value[1],
    shap_explanation[1],
    obs,
    feature_names=["feature 1", "feature 2"],
    matplotlib=True,
)
# Afficher les résultats LIME
lime_explanation.show_in_notebook()

divide by zero encountered in scalar divide

Commentaire :

Nous avons trouvé des explications cohérentes en comparant LIME et SHAP sur des données non linéaires générées par make_moons. Bien qu'ils n'affichent pas la même valeur, comme expliqué précédemment, ils attribuent le même poids aux caractéristiques pour chaque explication. Par exemple, pour la dernière cellule, c'est la feature 2 qui contribue à classer la donnée dans la classe 1 et la feature 1 qui contribue à classer la donnée vers la classe 0.

La cohérence des explications fournies par LIME et SHAP suggère que, pour les données non linéaires générées par make_moons, les deux méthodes sont capables de capturer et de refléter avec précision les modèles sous-jacents utilisés par le modèle pour ses prédictions. Cela pourrait être dû au fait que la frontière de décision, bien que non linéaire, est relativement simple et cohérente dans sa forme — des caractéristiques que LIME et SHAP peuvent efficacement saisir. Cela, ajouté au fait que nous n'avons que deux caractéristiques, peut simplifier l'interprétation et la comparaison des contributions de chaque caractéristique, rendant ainsi les explications de LIME et SHAP plus directement comparables. Si le modèle est bien ajusté et représente fidèlement la complexité sous-jacente des données, il est probable que les deux méthodes révèlent des tendances similaires dans la manière dont les caractéristiques influencent les prédictions, ce qui renforce la confiance dans la fiabilité et la validité des interprétations fournies.

Note : Ici, on peut voir que ce ne sont plus le log odds ratio qui est affiché mais bien la probabilité brute. Cela est du à comment l'Explainer SHAP (TreeExplainer) est construit... et comment SHAP décide de faire les choses. Par défaut, model_output="raw", et ici, avec notre random forest, SHAP décide d'utiliser la probabilité (alors qu'avec XGBoost, il utiliserait le log odds (des questions de choix d'implémentation..?)).

Générer d'autres explications globales

On travaille sur les données Adult (tâche de classification pour prédire si une personne gagne plus de 50 000 $ dans les années 1990).

Pas trop sûr de ce qui est demandé ici, alors j'ai choisi de m'amuser en générant plein de graphiques trop stylés de SHAP.

X, y = shap.datasets.adult()
X.head()

	Age	Workclass	Education-Num	Marital Status	Occupation	Relationship	Race	Sex	Capital Gain	Capital Loss	Hours per week	Country
0	39.0	7	13.0	4	1	0	4	1	2174.0	0.0	40.0	39
1	50.0	6	13.0	2	4	4	4	1	0.0	0.0	13.0	39
2	38.0	4	9.0	0	6	0	4	1	0.0	0.0	40.0	39
3	53.0	4	7.0	2	6	4	2	1	0.0	0.0	40.0	39
4	28.0	4	13.0	2	10	5	2	0	0.0	0.0	40.0	5

X_display, y_display = shap.datasets.adult(display=True)
X_display.head()

	Age	Workclass	Education-Num	Marital Status	Occupation	Relationship	Race	Sex	Capital Gain	Capital Loss	Hours per week	Country
0	39.0	State-gov	13.0	Never-married	Adm-clerical	Not-in-family	White	Male	2174.0	0.0	40.0	United-States
1	50.0	Self-emp-not-inc	13.0	Married-civ-spouse	Exec-managerial	Husband	White	Male	0.0	0.0	13.0	United-States
2	38.0	Private	9.0	Divorced	Handlers-cleaners	Not-in-family	White	Male	0.0	0.0	40.0	United-States
3	53.0	Private	7.0	Married-civ-spouse	Handlers-cleaners	Husband	Black	Male	0.0	0.0	40.0	United-States
4	28.0	Private	13.0	Married-civ-spouse	Prof-specialty	Wife	Black	Female	0.0	0.0	40.0	Cuba

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

tree = RandomForestClassifier(random_state=42)
tree.fit(X_train, y_train)
y_pred = tree.predict(X_test)
print("=" * 50)
print(f"Accuracy = {(y_pred == y_test).sum() / len(y_test):.4f}")
print("=" * 50)

==================================================
Accuracy = 0.8523
==================================================

explainer = shap.TreeExplainer(tree)
features = X_test[:20]
features_display = X_display[:20]
shap_values = explainer.shap_values(features)

shap.summary_plot(shap_values, features_display)

Le lien de parenté (Relationship) est étonnamment la caractéristique la plus importante.

expected_value = explainer.expected_value[1]
print(f"Explainer expected value : {expected_value}")

shap_values = shap_values[1]

Explainer expected value : 0.24102464680589678

shap.decision_plot(expected_value, shap_values, features_display)

misclassified = y_pred != y_test
misclassified = misclassified[:20]
print(f"Nombre d'exemples mal classifiés : {misclassified.sum()}")
shap.decision_plot(
    expected_value, shap_values, features_display, highlight=misclassified
)

Nombre d'exemples mal classifiés : 1

for (_, obs), shap_val in zip(
    features_display[misclassified].iterrows(), shap_values[misclassified]
):
    shap.decision_plot(
        expected_value,
        shap_val,
        obs,
        highlight=0,
    )

Commentaire decision plots :

On affiche un decision plot pour les 20 premières observations du jeu de test (moins lourd à calculer, et sinon ça fait beaucoup de lignes à afficher sur notre graphique, ce qui le rendrait illisible).

• L'axe des abscisses représente la sortie du modèle, qui est dans ce cas est la probabilité. Le graphique est centré sur la base value (expected_value), ici de 0.241.
• L'axe des ordonnées énumère les caractéristiques du modèle.

Toutes les valeurs SHAP sont relatives à la valeur attendue du modèle, comme les effets d'un modèle linéaire sont relatifs à l'ordonnée à l'origine. Par défaut, les caractéristiques sont classées par ordre décroissant d'importance. L'importance est calculée sur les observations représentées. Elle est généralement différente de l'ordre d'importance de l'ensemble des données.

La prédiction de chaque observation est représentée par une ligne colorée. En haut du graphique, chaque ligne touche l'axe des $x$ au niveau de la valeur prédite de l'observation correspondante. Cette valeur détermine la couleur de la ligne. En se déplaçant du bas du graphique vers le haut, les valeurs SHAP de chaque caractéristique sont ajoutées à la valeur de base du modèle. Cela montre comment chaque caractéristique contribue à la prédiction globale.

Ce graphique est notamment particulièrement informatif lorsque l'on cherche à comprendre pourquoi un exemple est mal classifié, ce que les deux derniers graphiques font.

Les graphiques d'après ne marchent pas avec un arbre de décision... :( Chaque explication marche donc plus ou moins bien en fonction du modèle utilisé !

On va donc utiliser un modèle de gradient boosting XGBoost... (en plus les valeurs de Shapley prennent beaucoup moins de temps à calculer !)

import xgboost

# train XGBoost model
model = xgboost.XGBClassifier(n_estimators=100, max_depth=2)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
print("=" * 50)
print(f"Accuracy = {(y_pred == y_test).sum() / len(y_test):.4f}")
print("=" * 50)

==================================================
Accuracy = 0.8746
==================================================

# compute SHAP values
explainer = shap.Explainer(model, X_train)
shap_values = explainer(X_test[:1000])

[16:36:19] WARNING: /workspace/src/c_api/c_api.cc:1240: Saving into deprecated binary model format, please consider using `json` or `ubj`. Model format will default to JSON in XGBoost 2.2 if not specified.

shap.plots.heatmap(shap_values)

Commentaire heatmap :

J'aime beaucoup ce graphique, il est très joli...

La sortie du modèle est affichée au-dessus de la matrice de la heatmap (centrée autour de la base value), et l'importance globale (valeur absolue) de chaque entrée du modèle est affichée sous la forme d'un diagramme à barres sur le côté droit du diagramme.

L'atout majeur de ce graphique réside dans sa capacité à illustrer la distribution de l'influence d'une caractéristique sur un sous-ensemble spécifique du jeu de données, offrant ainsi une perspective plus globale. Par exemple, il est notable que bien que les caractéristiques Capital Gain et Capital Loss influencent rarement le modèle, leur impact est particulièrement marqué lorsqu'elles sont effectivement présentes.

shap.plots.heatmap(shap_values, instance_order=shap_values.sum(1))

Par défaut, les exemples sont regroupés et ordonnés en fonction de la similarité de leurs explications. Il est possible de les trier en fonction de la somme de leurs valeurs de SHAP pour obtenir une autre perspective.

shap.plots.scatter(shap_values[:, "Age"])

Commentaire scatter :

Le graphique ci-dessus est un "dependence scatter plot". Il montre l'effet d'une caractéristique unique sur les prédictions faites par le modèle. Dans cet exemple, le logarithme des chances de gagner plus de 50 000 dollars augmente considérablement entre 20 et 40 ans.

• Chaque point représente une prédiction unique de l'ensemble de données.
• L'axe des $x$ est la valeur de la caractéristique (de la matrice $X$).
• L'axe des $y$ est la valeur SHAP pour cette caractéristique, qui représente dans quelle mesure la connaissance de la valeur de cette caractéristique modifie la sortie du modèle pour la prédiction de cet échantillon. Pour ce modèle, les unités sont les logarithmes des chances de gagner plus de 50 000 dollars par an.
• La zone gris clair au bas du graphique est un histogramme montrant la distribution des valeurs des données.

La dispersion verticale dans le graphique ci-dessus montre que la même valeur pour la caractéristique "âge" peut avoir un impact différent sur les résultats du modèle pour différentes personnes. Cela signifie qu'il existe des effets d'interaction non linéaires dans le modèle entre l'âge et d'autres caractéristiques.

shap_values.display_data = X_display.values
shap.plots.scatter(shap_values[:, "Age"], color=shap_values)

Commentaire scatter :

Pour montrer quelle caractéristique peut être à l'origine de ces effets d'interaction, nous pouvons colorer notre diagramme de dispersion de la dépendance à l'égard de l'âge par une autre caractéristique. Si un effet d'interaction est présent entre cette autre caractéristique et la caractéristique que nous traçons, il apparaîtra sous la forme d'un motif de coloration vertical distinct. Dans l'exemple ci-dessus, les jeunes de 20 ans ayant un niveau d'éducation élevé sont moins susceptibles de gagner plus de 50 000 dollars que les jeunes de 20 ans ayant un faible niveau d'éducation. Cela suggère un effet d'interaction entre le nombre d'années d'études et l'âge.

shap.plots.violin(shap_values)

Commentaire violon :

Le graphique du violon offre une représentation compacte de la distribution et de la variabilité des valeurs SHAP pour chaque caractéristique. C'est une autre façon d'afficher un summary plot ou un beeswarm. Les diagrammes en forme de violon sont empilés en fonction de l'importance de la caractéristique particulière sur la sortie du modèle (somme des valeurs absolues des valeurs SHAP par caractéristique).

Les diagrammes en forme de violon utilisent des figures en forme de violon pour afficher la distribution et la densité des valeurs SHAP pour leur caractéristique respective. Les violons peuvent donc donner un aperçu de l'étendue, de la variabilité, de l'asymétrie, de la symétrie et de la multimodalité de la distribution des valeurs SHAP pour un élément spécifique.

Le graphique récapitulatif des violons permet de comparer l'importance des caractéristiques. Des violons plus larges indiquent une densité plus élevée et des valeurs plus fréquentes, ce qui donne une idée de l'importance relative de chaque caractéristique par rapport aux résultats du modèle.

4 - Données textuelles [facultatif]

Apprendre un modèle de classification random forest sur les données fetch_20newsgroups puis utiliser LIME et SHAP pour obtenir des explications sur la classification d'exemples.

À l'aide de la documentation du dataset fetch_20newsgroups, apprendre un modèle random forest de classification bi-classe de texte.

from nltk.corpus import stopwords
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics import (
    ConfusionMatrixDisplay,
    accuracy_score,
    classification_report,
)

categories = ["sci.space", "sci.crypt"]
newsgroups_train = datasets.fetch_20newsgroups(
    subset="train", categories=categories, remove=("headers", "footers", "quotes")
)
newsgroups_test = datasets.fetch_20newsgroups(
    subset="test", categories=categories, remove=("headers", "footers", "quotes")
)

On utilise des stopwords afin d'éliminer les mots courants (tels que is, this, etc.) qui risqueraient d'introduire du bruit dans l'analyse textuelle et afin d'améliorer notre classifieur.

vectorizer = TfidfVectorizer(lowercase=False, stop_words=stopwords.words("english"))
train_vectors = vectorizer.fit_transform(newsgroups_train.data)
test_vectors = vectorizer.transform(newsgroups_test.data)

rf = RandomForestClassifier(random_state=42)
rf.fit(train_vectors, newsgroups_train.target)

RandomForestClassifier(random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

y_pred = rf.predict(test_vectors)
accuracy = accuracy_score(newsgroups_test.target, y_pred)
print(f"Accuracy: {accuracy:.2f}")

report = classification_report(
    newsgroups_test.target, y_pred, target_names=newsgroups_test.target_names
)
print("Classification Report:\n", report)

Accuracy: 0.88
Classification Report:
               precision    recall  f1-score   support

   sci.crypt       0.94      0.80      0.87       396
   sci.space       0.83      0.95      0.88       394

    accuracy                           0.88       790
   macro avg       0.88      0.88      0.88       790
weighted avg       0.88      0.88      0.88       790

fig, ax = plt.subplots(figsize=(10, 5))
ConfusionMatrixDisplay.from_predictions(newsgroups_test.target, y_pred, ax=ax)
ax.xaxis.set_ticklabels(categories)
ax.yaxis.set_ticklabels(categories)

[Text(0, 0, 'sci.space'), Text(0, 1, 'sci.crypt')]

Utilisez LIME et SHAP pour obtenir des explications sur les prédictions obtenues sur des exemples.

LIME : facile, très bien documenté pour les données textes.

from lime.lime_text import LimeTextExplainer
from sklearn.pipeline import make_pipeline

c = make_pipeline(vectorizer, rf)
explainer_lime = LimeTextExplainer(class_names=newsgroups_test.target_names)

SHAP : mal documenté pour des données textes entrainés sur des classifieurs classiques...

feature_names = vectorizer.get_feature_names_out()
# Transformer les données vectorisés en DataFrame pour SHAP
train_df = pd.DataFrame(train_vectors.todense(), columns=feature_names)
test_df = pd.DataFrame(test_vectors.todense(), columns=feature_names)
train_df.shape, test_df.shape

((1188, 26515), (790, 26515))

explainer_shap = shap.TreeExplainer(rf, train_df, feature_names=feature_names)
shap_values = explainer_shap.shap_values(test_df)
shap_values_ = explainer_shap(test_df)

 99%|===================| 1564/1580 [01:19<00:00]        

Première remarque : il faut un peu plus d'1 minute (soit environ 2 minutes 30 pour calculer la cellule du dessus) pour calculer les valeurs de Shapley sur les données de test, qui fait écho à la lenteur de SHAP, à cause de son approche "globale". Après, sachant qu'il y a 26515 caractéristiques (par 1188 données de train et 790 données de test) (taille du vocabulaire composant notre Bag of Words), cela reste acceptable.

idx = 69
exp = explainer_lime.explain_instance(
    newsgroups_test.data[idx], c.predict_proba, num_features=10
)
print(f"True class: {newsgroups_test.target_names[newsgroups_test.target[idx]]}")
print(f"Predicted class: {newsgroups_test.target_names[y_pred[idx]]}")
exp.show_in_notebook()

True class: sci.crypt
Predicted class: sci.crypt

exp.as_list()

[('Clipper', -0.10849659288809789),
 ('keys', -0.09476407916807551),
 ('phone', -0.08768716528119361),
 ('number', -0.03757703838836422),
 ('serial', -0.03512185928618738),
 ('The', 0.03162132353993653),
 ('get', 0.026208618806981658),
 ('care', -0.0247021558381579),
 ('warrant', -0.022426825823142052),
 ('using', -0.012622610353977284)]

Commentaire LIME :

La visualisation LIME montre que les termes "Clipper", "keys", "phone", "number", et "serial" sont les plus influents pour classer un document spécifique comme relevant de la catégorie sci.crypt. Les mots "care", "warrant", et "using" ont également une contribution positive, bien que moins significative. Cependant, les mots "The" et "get" ont une influence négative, favorisant à tort la classification dans la catégorie sci.space. Malgré cette influence négative, elle est insuffisante pour renverser la tendance imposée par les termes clés positifs.

class_predite = y_pred[idx]
shap.plots.waterfall(shap_values_[idx, :, class_predite])

Commentaire SHAP :

Nous avons choisi d'utiliser le plot "waterfall" pour la visualisation. Il est similaire au "force plot" mais affiche les caractéristiques d'une manière cascadée. Cela rend le graphique plus lisible, compte tenu de la forme de nos caractéristiques. De plus, on remarque que c'est bien la probabilité d'obtenir la classe qui est affichée. Ici, $ P(\text{sci.crypt}) = 0.82 $, comme avec LIME (normal, vu que l'on travaille sur le même classifieur).

Note : Dans SHAP, on doit préciser, en plus de l'exemple, quelle classe on souhaite expliquer, c'est-à-dire la classe prédite par le classifieur ; ici 0 (sci.crypt). On peut, si on le souhaite, expliquer l'autre classe, mais, ce n'est généralement pas ce que l'on veut faire en explicabilité.

Ainsi, le graphique confirme que "keys" est le terme le plus déterminant pour la classification en sci.crypt, suivi de près par "Clipper" et "phone". Des termes tels que "number", "serial", et "warrant" jouent également un rôle important. Les termes "tap", "using", et "care" apparaissent comme moins déterminants mais sont tout de même des indicateurs positifs de la classe sci.crypt. Il est intéressant de noter que SHAP attribue une importance relative légèrement différente à certains mots par rapport à LIME, mais les deux méthodes s'accordent globalement sur les caractéristiques les plus influentes. SHAP révèle également l'influence combinée des autres caractéristiques, montrant que leur contribution globale est mineure par rapport aux termes clés identifiés.

Les deux visualisations mettent en lumière les mots-clés influençant le modèle, mais la cohérence entre les deux méthodes renforce la confiance dans l'interprétation des caractéristiques importantes. Les différences dans l'ordre et l'importance des caractéristiques reflètent les nuances des deux méthodes d'explication. D'un côté, LIME fournit une approximation locale linéaire de la prédiction, qui peut être plus interprétable mais moins précise dans le cas de modèles non linéaires. De l'autre, SHAP offre une perspective basée sur la théorie des jeux et considère les contributions de toutes les caractéristiques ensemble, ce qui peut donner une vue plus globale mais potentiellement plus complexe à interpréter.

idx = -1
exp = explainer_lime.explain_instance(
    newsgroups_test.data[idx], c.predict_proba, num_features=6
)
print(f"True class: {categories[newsgroups_test.target[idx]]}")
print(f"Predicted class: {newsgroups_test.target_names[y_pred[idx]]}")
exp.show_in_notebook()

True class: sci.space
Predicted class: sci.space

class_predite = y_pred[idx]
shap.plots.waterfall(shap_values_[idx, :, 1])

Ci-dessus, on peut faire le même commentaire que précédemment. On retrouve des similitudes entre les explications données par SHAP et LIME. Cependant, une remarque intéressante et à quel point les 26506 autres features jouent un role.

Par contre, contrairement à LIME, SHAP nous permet d'avoir une explication "globale" :

shap.summary_plot(shap_values, test_df, plot_type="bar", feature_names=feature_names)

Commentaire :

L'analyse du summary plot révèle une première observation intrigante : chaque terme semble avoir une contribution similaire pour les différentes classes. Par exemple, le terme "cost" montre une valeur absolue moyenne de SHAP identique pour les classes 1 et 0, soit 0.005. Cette découverte met une nouvelle fois en lumière une contrainte de ce graphique, lequel ne permet pas d'appréhender la direction de l'influence des termes.

Toutefois, en examinant de plus près les graphiques suivants pour chacune des classes, la raison de cette particularité se clarifie. Il apparaît que si un terme a une forte influence (i.e. une forte valeur) pour classer un texte dans une catégorie donnée, il possède un pouvoir de classification équivalent mais opposé pour l'autre catégorie. Cela est illustré par le terme "encryption", qui joue un rôle prépondérant pour attribuer un texte à la classe sci.crypt. En contrepartie, une valeur très basse pour ce terme pèse d'autant plus dans la balance pour orienter la classification vers l'autre classe, sci.space.

shap.summary_plot(shap_values[0], test_df, feature_names=feature_names)  # sci.space

shap.summary_plot(shap_values[1], test_df, feature_names=feature_names)  # sci.crypt